function [params,fitdata] = invffitfn(data,showplot,numfitdatasamples,init)

% function [params,fitdata] = invffitfn(data,showplot,numfitdatasamples,init)
%
% fits a/(f^b) (pink noise or 1/f function) to data. 
% data: frequency (column 1) and amplitude (column 2)
% numfitdatasamples: number of samples applied to the fitting function
% init: initial parameter values [a,b]
% showplot: flag to display fit in graph (in log-log coorinates)
%
% params params(1) = scaling factor, params(2) = exponent / slope in
% log-log coordinates
% fitdata: 2 column data output that can be plotted to show fit
%
% May 2017  JMG

defarg('showplot',1);
defarg('numfitdatasamples',1000);
defarg('init',[10000 1.5]);
sz = size(data);
if sz(2) ~= 2
	params = [];
    fitdata = [];
    disp('Data incomplete.');
    return
end

% y = a/(f^b)
errfn = inline('sum(((x(1)./(P1(:,1).^x(2))- P1(:,2)).^2))',1);
options = optimset('MaxIter',100000,'TolFun',0.0000001,'TolX',1e-7,'MaxFunEvals',20000);
params = fminsearch(errfn,init,options,data);

% plot
fitdata = linspace(min(data(:,1)),max(data(:,1)),1000)';
fitdata = [fitdata,params(1)./(fitdata.^params(2))];
if showplot,
    figure;
    loglog(data(:,1),data(:,2),'k*'); hold on;
    loglog(fitdata(:,1),fitdata(:,2),'b:');
    drawnow;
end

return
